function pde = Kelloggdata(varargin)
%%  KELLOGGDATA data of Kellogg Problem
%  u = r^{gamma} sin(mu \theta)
%  theta is the polar angle either in cylindrical or polar coords
%  Reference Z. Chen and S. Dai 2002 SISC
%  "On the efficiency of adaptive finite element methods for elliptic problems
%   with discontinuous coefficients"
%
% Copyright (C) Long Chen. See COPYRIGHT.txt for details.

rho = pi/4;

if isempty(varargin) || varargin{1} == 0.1
    gamma = 0.1;
    sigma = -14.92256510455152;
    R = 161.4476387975881;
elseif varargin{1} == 0.5
    gamma = 0.5;
    sigma = -2.3561944901923448;
    R = 5.8284271247461907;
elseif varargin{1} == 0.02
    gamma = 0.02;
    sigma = -77.754418176347386;
    R = 4052.1806954768103;
end


pde = struct('f',0,...
    'exactu',@exactu,...
    'g_D',@exactu,...
    'Du',@Du,...
    'd',@d);

    function K =  d(p)  % Diffusion constant
        idx = (p(:,1).*p(:,2) >0);
        K  = ones(size(p,1),1);
        K(idx) = R;
    end

    function u =  exactu(p) % exact solution
%         r = sqrt(sum(p.^2,2)); % this will be problematic in 3D
        r = sqrt(p(:,1).^2+p(:,2).^2); % good for both 2D and 3D
        theta = atan2(p(:,2),p(:,1));
        theta = (theta>= 0).*theta + (theta<0).*(theta+2*pi);
        mu = (theta>= 0 & theta<pi/2).*cos((pi/2-sigma)*gamma).*cos((theta-pi/2+rho)*gamma)...
            +(theta>= pi/2 & theta<pi).*cos(rho*gamma).*cos((theta-pi+sigma)*gamma)...
            +(theta>= pi & theta<1.5*pi).*cos(sigma*gamma).*cos((theta-pi-rho)*gamma)...
            +(theta>= 1.5*pi & theta<2*pi).*cos((pi/2-rho)*gamma).*cos((theta-1.5*pi-sigma)*gamma);
        u =  r.^gamma.*mu;
    end

    function uprime =  Du(p)
        % the gradient of exact solution
        theta = atan2(p(:,2),p(:,1));  % polar angle
        theta = (theta>= 0).*theta +(theta<0).*(theta+2*pi);
        t = 1+(p(:,2).^2)./(p(:,1).^2);
        % r = sqrt(sum(p.^2,2));
        r = sqrt(p(:,1).^2+p(:,2).^2); % good for both 2D and 3D
        rg = r.^gamma;
        
        ux1 = (p(:,1)> 0.0 & p(:,2)>= 0.0).*...
            (rg.*gamma./r.*cos((pi/2-sigma)*gamma)./r.*p(:,1)...
            .*cos((theta-pi/2+rho)*gamma)...
            +rg.*cos((pi/2-sigma)*gamma).*sin((theta-pi/2+rho)*gamma)...
            *gamma.*p(:,2)./(p(:,1).^2)./t);
        
        uy1 = (p(:,1)> 0.0 & p(:,2)>= 0.0).*...
            (rg*gamma./r.*cos((pi/2-sigma).*gamma)...
            .*cos((theta-pi/2+rho).*gamma)./r.*p(:,2)...
            -rg.*cos((pi/2-sigma).*gamma).*sin((theta-pi/2+rho)*gamma)...
            *gamma./p(:,1)./t);
        
        ux2 = (p(:,1)<= 0.0 & p(:,2)> 0.0).*...
            (rg*gamma./r.*cos(rho.*gamma)./r.*p(:,1)...
            .*cos((theta-pi+sigma).*gamma)...
            +rg.*cos(rho*gamma).*sin((theta-pi+sigma).*gamma)*gamma...
            .*p(:,2)./(p(:,1).^2)./t);
        
        uy2 = (p(:,1)<= 0.0 & p(:,2)> 0.0).*...
            (rg*gamma./r.*cos(rho.*gamma)./r.*p(:,2)...
            .*cos((theta-pi+sigma)*gamma)-rg.*cos(rho*gamma)...
            .*sin((theta-pi+sigma).*gamma)*gamma./p(:,1)./t);
        
        ux3 = (p(:,1)< 0.0 & p(:,2)<= 0.0).*...
            (rg*gamma./r.*cos(sigma.*gamma)./r.*p(:,1)...
            .*cos((theta-pi-rho).*gamma) ...
            +rg.*cos(sigma.*gamma).*sin((theta-pi-rho).*gamma).*gamma...
            .*p(:,2)./(p(:,1).^2)./t);
        
        uy3 = (p(:,1)< 0.0 & p(:,2)<= 0.0).*...
            (rg*gamma./r.*cos(sigma.*gamma)./r.*p(:,2)...
            .*cos((theta-pi-rho).*gamma) -rg.*cos(sigma*gamma)...
            .*sin((theta-pi-rho).*gamma)*gamma./p(:,1)./t);
        
        ux4 = (p(:,1)>= 0.0& p(:,2)< 0.0).*...
            (rg*gamma./r.*cos((pi/2-rho).*gamma)./r.*p(:,1)...
            .*cos((theta-3*pi/2-sigma)*gamma) ...
            +rg.*cos((pi/2-rho).*gamma).*sin((theta-3*pi/2-sigma)*gamma)...
            *gamma.*p(:,2)./(p(:,1).^2)./t);
        
        uy4 = (p(:,1)>= 0.0 & p(:,2)< 0.0).*...
            (rg*gamma./r.*cos((pi/2-rho).*gamma)./r.*p(:,2)...
            .*cos((theta-3*pi/2-sigma)*gamma)...
            -rg.*cos((pi/2-rho)*gamma).*sin((theta-3*pi/2-sigma)*gamma)...
            *gamma./p(:,1)./t);
        
        uprime(:,1) = ux1+ux2+ux3+ux4;
        uprime(:,2) = uy1+uy2+uy3+uy4;
        if size(p,2) == 3
            uprime(:,3) = zeros(size(p,1),1);
        end
    end
end